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I.  Introduction 

Issues  concerning  LiH  thermal  decomposition  and  LiH  and  Li  reactions  with  air,  discussed  in  a  general  form  in  [1], 
have  now  to  be  investigated  within  the  scramjet  combustion  chamber.  Here,  thanks  to  the  compression  operated  by  the 
oblique  shock  system,  the  supersonic  flow  posses  high  temperature  and  high  dynamic  pressure;  entering  into  the 
combustion  chamber  it  flows  over  the  solid  grain  and,  as  an  effect  of  the  shear  stresses  on  the  surface,  a  turbulent 
compressible  boundary  layer  grows  along  the  chamber  while  the  grain  surface  regresses;  the  mass  produced  by  the  fuel 
consumption  is  injected  into  the  boundary  layer  and  increases  its  thickness.  Thus  inside  the  boundary  layer  and/or  into 
the  core  flow  the  fuel  must  react  with  the  hot  air  stream  in  a  few  milliseconds  increasing  the  flow  temperature  and 
producing  the  thermal  fluxes  necessaries  to  gasify  the  solid  grain.  For  these  reasons  critical  issues  are:  species  diffusion 
into  the  hot  stream,  and  their  reaction  rates  with  air;  these  issues  in  general  make  gaseous  hydrogen  preferable  as  fuel 
for  scramjet  application  (because  its  chemical  and  diffusion  times  are  shorter  than  those  of  other  fuels  when  compared 
to  convective  time).  From  this  viewpoint  the  particular  behaviour  of  LiH,  that  can  release  gaseous  hydrogen  by  thermal 
decomposition,  together  with  light  and  highly  reactive  lithium,  offers  an  ideal  solution  to  both  hydrogen  storage  and 
high  energy  density  problems. 

Aim  of  this  report  is  to  present  the  physical  model  describing  LiH  “combustion”,  and  to  discuss  by  means  of  some 
approximated  models  the  main  topics  associated  with  thermal  fluxes,  diffusion  rate,  lithium  vaporization  and 
combustion.  To  give  a  rough  numerical  estimate  of  the  characteristic  times  involved  a  “flight”  configuration  was 
investigated  similar  to  that  proposed  by  Jarymovycz,  Yang  and  Kuo  and  described  in  [1]. 


I.l  The  Physical  Model 

As  a  consequence  of  lithium  hydride  and  lithium  thermochemical  properties  examined  in  [2]  and  of  the  flow  conditions 
into  the  combustion  chamber  summarized  above,  hereinafter  we  will  assume  (justifying  each  assumption  in  the 
following  sections)  the  grain  “combustion”  to  be  described  by  the  following  steps: 

Due  to  the  thermal  fluxes  produced  by  combustion  in  the  gas  phase  the  LiH  grain  surface  liquefies;  dealing 
with  a  phase  change  we  can  assume  the  interface  between  solid  and  liquid  LiH  (or  “wall”)  to  be  always  at  the 
LiH  liquefaction  temperature  (  T^q  =  950  K  at  1  atmosphere  ). 

Thanks  to  liquid  LiH  high  thermal  diffusivity,  the  small  LiH  droplets  formed  at  the  wall  [Litfire],  heated  by  the 
thermal  fluxes  reach  the  LiH  decomposition  temperature  (Tjec  =  1100  K);  LiH  decomposes  fast  when 
compared  with  the  combustor  convective  time  and  before  reacting  with  oxygen,  its  decomposition  time  being 
of  the  order  of  10'^  s;"*^’  in  fact  the  LiH  reaction  rates  with  oxygen  is  lower  than  its  decomposition  rate  (at 
least  by  three  orders  of  magnitude)  and  the  oxygen  concentration  near  the  surface  is  too  small  to  allow  any 
kind  of  oxidation.^”*’  As  stated  in  [2],  LiH  does  not  react  with  air  nitrogen  in  this  range  of  temperatures. 

As  an  effect  of  the  thermal  decomposition,  gaseous  hydrogen  leaves  the  drop  surface  bubbling  vigorously  and 
contributing  to  its  fragmentation  into  smaller  liquid  lithium  droplets.  Hydrogen  diffuses  through  the  boundary 
layer  and  mixed  with  air  oxygen  by  the  highly  turbulent  stream,  burns  and  releases  heat.  Some  lithium  particles 
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may  be  trapped  into  the  hydrogen  bubbles  leaving  the  surface,  burning  with  oxygen  and  increasing  the  heat 
released  by  hydrogen  combustion. 

Liquid  Li  droplets  are  dragged  away  from  the  surface  (as  an  effect  of  the  intense  shear  stresses  due  to  the  high 
speed  flow)  into  the  turbulent  boundary  layer. 

The  light  liquid  lithium  droplets  thus  are  mixed  with  air  and  hydrogen  combustion  products  present  into  the 
boundary  layer;  passing  through  the  hot  gas  in  the  core  flow  they  are  vaporized  by  the  thermal  fluxes,  yielding 
highly  reactive  gaseous  lithium  that,  transported  by  the  turbulent  stream,  burns  with  air  oxygen  or  with  the 
radicals  produced  by  the  hydrogen  combustion  (see  [2]).  When  allowed  by  the  local  oxygen  concentration, 
lithium  combustion  occurs  in  proximity  of  the  droplet  surface,  enhancing  liquid  lithium  vaporization.  Again,  as 
stated  in  [2],  liquid  lithium  doesn’t  react  with  air  nitrogen  in  this  range  of  temperatures. 

Where  the  combustion  is  well  developed  the  water  vapour  concentration  is  so  high  that  it  can  react  with  the 
liquid  LiH  near  the  surface;  however  the  net  effect  of  this  reaction  is  the  production  of  LiOH,  which  at  these 
temperatures  decomposes  yielding  liquid  lithium  and  OH. 


1.2  Inside  the  boundary  layer 

The  first  step  to  assess  the  validity  of  this  physical  model  is  the  evaluation  of  thermal  fluxes  and  shear  stresses  acting  on 
the  grain  surface.  At  the  combustion  chamber  entrance,  thanks  to  the  compression  operated  by  the  oblique  shock 
system,  the  air  flow  is  at  high  temperature;  thus  we  are  in  presence  of  high  stagnation  enthalpy  over  the  grain  that  can 
heat  the  surface  beyond  the  LiH  liquefaction  point.  We  therefore  begin  by  investigating  the  heat  transfer  rate  magnitude 
assuming  non  reacting  flow  conditions.  The  idea  is  that  if  the  convective  thermal  fluxes  due  solely  to  flow  conditions 
can  promote  the  grain  surface  liquefaction,  the  heat  released  by  the  combustion  surely  will  be  able  to  feed  a  self- 
sustained  combustion  process. 

The  simplest  preliminary  estimate  of  the  convective  heat  transfer  between  a  wall  and  a  fluid  with  variable  properties  is 
provided  by:®^ 


q*  =St-p,-V^-(h^„-h„) 


(1.1) 


where  and  Vg  are,  respectively,  the  flow  density  and  (relative)  velocity  at  the  outer  edge  of  the  boundary  layer 

and  St  is  the  Stanton  number,  a  dimensionless  heat  transfer  coefficient  (representing  a  nondimensional  thermal 
conductance).  Finally  h^*  (adiabatic  wall  enthalpy  or  recovery  enthalpy)  is  the  stagnation  enthalpy  that  the  wall- 
bounding  streamline  would  attain  if  the  wall  were  perfectly  insulated  (the  corresponding  temperature  is  known  as  the 
adiabatic  wall  temperature  or  recovery  temperature).  In  fact,  if  conductive  heat  transfer  within  the  fluid  is  neglected,  a 
first  law  balance  shows  that  h^w  is  equal  to  the  stagnation  enthalpy  of  the  flow:  the  free  stream  kinetic  energy  is 
completely  converted  by  the  viscosity  into  internal  (thermal)  energy.  However,  because  of  the  small  but  finite  thermal 
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conductivity  of  the  fluid,  thermal  energy  is  conducted  away  from  the  stagnation  point  to  the  adjacent  gas,  thus  lowering 
the  value  of  h;,*;  to  account  for  this  phenomenon  it  is  customary  to  introduce  a  recovery  factor  r  (r  <1): 

haw=r~  +  he(TJ  (1.2) 

The  term  h„  (eqn  1. 1)  is  the  static  enthalpy  at  the  surface,  corresponding  to  the  design  value  for  cooled  surfaces  or  to 
that  of  the  propellant  at  the  phase  change  temperature.  Depending  on  its  value  (if  higher  or  lower  than  h^w)  it  is  possible 
to  determine  the  thermal  flux  direction  (figure  1). 


Fig.  1:  Convective  heat  fluxes  on  a  surface  where  V  =  0 


To  determine  the  adiabatic  wall  enthalpy,  the  recovery  factor  and  the  flow  properties  into  the  boundary  layer,  our 
analysis  must  account  for  the  flow  compressibility;  in  fact  the  presence  of  high  velocities  in  the  boundary  layer  gives 
rise  to  such  large  temperature  differences  that  it  is  necessary  to  account  for  the  changes  on  fluid  properties  together 
with  that  on  volume.  Furthermore  the  heat  transfer  plays  an  important  part  in  the  boundary  layer  behaviour  leading  to 
the  appearance  of  strong  interaction  between  the  velocity  field  and  the  temperature  field.  As  known  in  literature,  it  is 
possible  to  describe  the  compressible  boundary  layer  by  means  of  the  same  relations  holding  for  incompressible  flow 
(for  both  laminar  and  turbulent  cases)  on  condition  that  the  value  of  density  and  viscosity  are  taken  at  a  suitably  chosen 
reference  temperature  T*.®^’  There  exist  in  literature  various  methods  to  determine  T*;  for  our  rough  analysis  the 
reference  temperature  and  the  heat  flux  can  be  evaluated  for  a  turbulent  flow  on  a  flat  plate  by  the  Eckert  reference 
enthalpy  method,  assuming  zero-pressure  gradient,  non  reacting  surface  and  constant  wall  temperature.^^ 


1.2.1  The  Eckert  reference  enthalpy  method** 

This  method  is  based  on  the  assumption  that  a  variable  properties  boundary  layer  can  be  cast  in  a  form  similar  to  that 
of  a  constant  properties  boundary  layer  when  computed  for  a  proper  “reference”  enthalpy  (and  temperature); 
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hereinafter  the  “reference”  quantities  used  in  estimating  will  be  designated  with  the  superscript  First  the 
reference  enthalpy  is  found  by  solving  simultaneously  the  following  equations: 


h  = 


h  +h„ 


■  +  0.22-r 


V. 


and 


(1.3) 


(1.4) 


This  solution  is  iterative  due  to  the  interdependence  between  h*  (or  T*)  and  the  air  properties.  After  calculating  the 
reference  h*  and  the  corresponding  reference  temperature  T*,  the  Stanton  number  for  zero  pressure  gradient,  constant 
wall  temperature  (and  for  a  flat  plate  turbulent  flow)  is: 


with: 


0.0296 

(Pry  (Re:)5 


Re  = 


P*  VeX 

* 

P 


(1.5) 


(1.6) 


* 

P 


Pe 

RT* 


and 


Pr*  = 


In  this  simplified  analysis  the  convective  heat  transfer  at  the  interface  is  thus: 


(1.7) 


(1.8) 


:St*  P*V,  (h,.-hj 


(1.9) 


At  the  same  temperature  T*  we  can  estimate  other  boundary  layer  properties  such  as  its  thickness  8*(x): 


5*(x)  0.370 


(1. 10) 
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and  the  turbulent  skin  friction  coefficient  (Reynolds  analogy): 


* 


0.370 


=  2St*(Pr*)3 


(I.ll) 


1.2.2  Convective  heat  transfer  and  shear  stresses 

As  stated  in  the  introduction,  the  order  of  magnitude  of  convective  heat  flux  and  shear  stress  have  been  examined  using 
as  reference  case  a  “flight”  configuration  similar  to  that  proposed  by  Jarymovycz,  Yang  and  K.  Kuo.^"*  In  our  case  we 
have  assumed  the  flow  entering  into  a  two-dimensional  planar  combustion  chamber  with  a  Mach  number  M  =  2,  static 
pressure  Pe  =  1  atm  and  static  temperature  T  =  1400  K,  corresponding  to  a  free  stream  velocity  Vg  =  1451  m/s.  The 
chamber  length  was  assumed  to  be  L  =  1  m,  its  height  about  20  cm  while  the  LiH  grain  surface  spans  50  cm  along  the 
lower  boundary  of  the  chamber,  thus  allowing  space  for  chemical  reactions  downstream  of  the  fuel  grain. 


50  cm 

1 

Fig.  2  :  Schematic  view  of  the  combustion  chamber  used 


The  Air  and  LiH  properties  (enthalpies,  constant  pressure  heat  ratio,  viscosity  and  thermal  conductivity)  were  calculate 
using  the  NASA  CEA600  database  and  extrapolating  data  from  the  Dry  Air  Properties  Table.® 

As  a  first  step  the  reference  enthalpy  h*  (and  thus  the  corresponding  T*)  was  calculated  using  the  Eckert  method;  the 
first  tentative  Ti*,  useful  to  evaluate  at  the  first  iteration  the  air  properties  in  (1.4),  was  calculated  using  the  Dorrance 
formula: 
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fry  y 

1  +0.032  M"  +0.58 

e 

iTe  ■  )_ 

(1.12) 


where  Te  and  Me  are  respectively  the  static  temperature  and  the  Mach  number  of  the  flow  at  the  edge  of  the  boundary 
layer  while  T*  is  the  LiH  surface  temperature  (  Tw  =  950  K)  assumed  constant,  since  it  is  the  phase  change  temperature. 
Thus  using  (1.9)  and  (1.1 1)  and  remembering  that: 


T 


W 


(1.13) 


the  heat  fluxes  and  the  shear  stress  along  the  grain  surface  were  calculated;  results  are  shown  in  figure  3. 
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Fig.  3:  Heat  fluxes  and  shear  stresses 


The  convective  heat  fluxes  and  shear  stresses  due  only  to  the  energy  recovery  from  the  flow  entering  and  passing  on  the 
surface  are  high  and  follow  the  boundary  layer  profile.  To  assess  the  effect  of  these  thermal  fluxes  on  the  surface,  a 
rough  estimate  of  the  lithium  hydride  liquefaction  rate  was  done.  We  assume  that  the  heat  absorbed  locally  by  the 
interface  is  the  sum  of  the  latent  heat  of  liquefaction  (see  table  1,  where  a  comparison  is  made  with  other  fuels  used  in 
literature  for  studies  of  solid  fuelled  scramjets  (see  [1]))  and  the  heat  required  to  increase  the  liquid  LiH  temperature 
from  950  K  up  to  1 100  K  (AH  =  875  kJ/kg). 


Fuel 

Average 

molecular 

formula 

[-] 

Mass 

density 

[kg/m'] 

Molar  mass 
[kg/kmol] 

Heat  of 
formation 
[kJ/mol] 

Melting 

point 

[K] 

Thermal 

conductivity 

[W/(m-K)] 

Heat  of 
gasification 
[kJ/kg] 

Plexiglass  or 
PMMA 

(C5H«02)„ 

1180 

100 

-430,5 

527.6 

0.17-0.19 

1300-2700 

PolyEthylene 

(C2H4)„ 

910-965 

28 

-53,8 

400 

0.40-0.46 

3000-5600 
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(PE) 

Polystyrene 

(PS) 

(CsHs)„ 

1050 

104 

+  18,4 

510 

0.1-0.13 

2700 

HTPB 

(CloHi5.5380o,073)n 

930 

138 

-51,8 

N/A 

0.217 

1800 

Lithium 

Hydride 

LiH 

800 

7.9 

-90.63 

950 

7.1 

2625* 

*  latent  heat  of  liquefaction 

Table  1:  Comparison  among  fuels  for  scrj  applications 


In  fact,  since  liquid  LiH  thermal  diffusivity  is  D*  ~  10’^  m^/s,  the  characteristic  heating  time  of  a  l)a.m  deep  liquid  layer 
is  about  10’®  s;  thus  we  can  assume  the  liquefaction  and  heating  processes  to  happen  by  means  of  a  unique  continuous 
step  with  a  total  local  heat  required  of  about  3400  kJ/kg.  So,  for  a  thermal  flux  q*  ~  1  MW/m^ ,  corresponding  to  q*  =  1 
J/s*mm^,  since  the  liquid  LiH  density  pi  =5x10'^  kg/mm^  the  local  flowrate  is  m^,  ~  0.6  mmVs  corresponding  to  a 
local  regression  rate  of  the  surface  of  about  0.6  mm/s;  since  the  convective  time  is  of  order  10'^  s,  the  thickness  of  the 
liquid  layer  formed  during  one  millisecond  on  the  surface  is  hi  =:  0.6pm.  This  result  will  be  useful  in  the  following 
section  to  assess  (together  with  other  considerations)  the  size  of  the  droplets  leaving  the  surface. 

Notice  that  the  surface  regression  rate  (corresponding  to  the  convective  heat  flux  just  calculated)  is  of  the  order  0.6 
mm/s;  this  value  is  slightly  lower  than  typical  regression  velocities  of  fuels  commonly  used  in  space  application. 
However  in  this  “conservative”  analysis  the  heat  flux  contribute  due  to  gas  phase  chemical  reactions  was  neglected;  this 
last  will  enhance  strongly  the  combustion  performance  by  increasing  the  regression  rate  and  thus  the  amount  of  gaseous 
fuel  produced.  An  example  of  the  thermal  power  available  if  all  the  fuel  mass  flow  rate  calculated  (  m^  =  0.6  mmVs) 
reacts  in  the  gas  phase  is  given  in  table  2  : 


Reoctnnts 

Pioducts 

If 

Theiinctl  Powei 
[W] 

Li  +  02 

Li20 

242 

10.9 

Li  +  02H 

Li20  +  H 

435 

19.6 

Li  +  02H 

LiOH  +  02 

600 

27 

Li  +  OH 

LiOH 

633 

28.5 

H2  +02 

H20 

253 

5.7 

Table  2:  Thermal  powers  corresponding  to  my 

where  the  thermal  power  obtainable  by  the  lithium  in  my  (assumed  to  react  completely  following  one  among  the 
presented  reactions)  is  compared  to  that  of  the  reaction  among  the  hydrogen  mass  contained  in  my  and  the  air  oxygen. 


1.2.3  Species  diffusion 

As  a  result  of  the  thermal  decomposition  gaseous  hydrogen  is  released  from  near  the  surface  and  is  transported  by  the 
turbulent  stream;  in  the  same  way  liquid  lithium  drops  exposed  to  the  high  temperature  flow  vaporize  releasing  gaseous 
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lithium.  In  the  following  sections  the  contribute  of  both  molecular  and  turbulent  mass  transport  will  be  valued  and 
compared  to  assess  the  ability  of  these  gaseous  species  to  diffuse  through  the  non  reacting  boundary  layer  up  to  the  core 
flow. 


1.2.3. 1  Molecular  diffusion 

As  a  first  step  the  convective  characteristic  time  t^  was  compared  with  both  hydrogen  and  gaseous  lithium  diffusive 
times  defined  as: 


t 


d 


(1.14) 


where  5  is  the  boundary  layer  thickness  ad  Dy  is  their  binary  diffusivity  into  the  air.  Data  on  hydrogen  diffusivity  are 
available  in  [Svelha]^**  while  those  concerning  lithium  were  calculated  (assuming  the  flow  to  be  nonreacting,  f.i.,  the  air 
to  be  composed  of  gaseous  oxygen  only)  using  the  following  relation:’*’’^ 


(1.15) 


where  ke  =  1.38  x  10'^^  [J/K]  is  the  Boltzmann  constant  and  r;  is  the  Li  atomic  radius  or  two  times  the  Van  der  Waals 
radius  for  H2  and  N2  molecules  while  mi  is  the  mass  of  gaseous  species.  To  verify  the  correctness  of  the  values 
obtained,  equation  1.15  was  used  also  to  calculate  the  hydrogen  diffusivity  which  was  successively  compared  with  [70]. 
Molecular  masses  and  radii  used  are  reported  in  table  3: 


Li 

H2 

02 

ma.ss  [kg*10^‘] 

1.146 

0.332 

4.65 

1X1(1111.9  [pin] 

182 

240 

304 

Table  3 

Alternatively  Kuo  suggests  a  method  based  on  the  bifurcation  approximation  of  Bird,’^  used  to  evaluate  the  molecular 
diffusion  flux  in  the  Reynolds-averaged  Navier-Stokes  equations;  in  this  method  the  diffusion  coefficient  Dij  is  replaced 
by 


D  (t  ,  p) 

g.  gj 


(1.16) 
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where  D  is  a  reference  diffusion  coefficient  given  by: 


1.719  X  10^^ 

P 


(1.17) 


where  D  is  in  cm7sec,  T  is  in  Kelvin  and  p  is  in  atmospheres.  The  gi  coefficients  are  referred  to  as  diffusion  factor  for 
species  i  and  can  be  calculated  from: 


g.  = 


v26.7 


(1.18) 


For  most  gas  systems,  these  correlation  are  within  5%  of  more  exact  values  for  temperatures  up  to  3000  K.  Both 
methods  give  results  which  are  in  good  agreement  (for  the  hydrogen  case)  with  [70];  hydrogen  and  lithium  diffusivities 
calculated  using  equation  ( 1. 16)  are  plotted  in  figure  4. 


Fig.  4:  Diffusivities  as  a  function  of  temperature 


Thus  values  corresponding  to  the  non  reacting  flow  temperature  (T  =  1400  K;  this  value  was  used  as  conservative  lower 
limits)  were  used  to  calculate  the  characteristic  diffusion  times;  results  are  shown  in  figure  5 
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Fig.  5:  Molecular  diffusion  times  and  boundary  layer  thickness  along  the  chamber  axis 


As  stated  before  this  is  the  minimum  time  that  both  hydrogen  and  lithium  take  to  diffuse  through  a  non  reacting 
boundary  layer  along  the  combustion  chamber  axis;  hydrogen  and  lithium  take  a  time  lower  than  convective  time  only 
in  the  first  part  of  the  grain  (about  1/5  of  the  grain  length);  however  the  species  diffusion  is  enhanced  by  the  increase  in 
temperature  due  to  the  combustion;  furthermore  while  hydrogen  diffuses  from  the  surface,  lithium  drops  vaporize  (see 
next  sections)  on  the  way;  gaseous  lithium  thus  can  diffuse  while  the  vaporizing  lithium  particle  is  near  the  upper  border 
of  the  boundary  layer,  so  reaching  immediately  the  core  flow. 


I.2.3.2  Turbulent  mass  transport 

A  simplest  estimate  of  the  turbulent  mass  transport  due  to  the  turbulence  was  made  using  the  concept  of  eddy  mass 
diffusivity  Dj.  As  well  known  in  literature,  when  dealing  with  the  time  (Reynolds)  averaged  Navier-Stokes  equations 
the  eddy  viscosity  Vj,  the  turbulent  thermal  diffusivity  aj  and  the  eddy  mass  diffusivity  Dj  are  used  to  express  the 
correlations  among  fluctuating  velocities,  enthalpy  and  species  in  terms  of  mean  flow  quantities  so  achieving  the 
equation  system  closure.  In  particular  the  eddy  mass  diffusivity  is  connected  to  the  eddy  viscosity  Vj  by  means  of  the 
turbulent  Schmidt  number  Sct: 


Dj. 


(1.19) 


In  most  of  models  developed  to  describe  the  evolution  of  turbulent  boundary  layers  the  turbulent  Schmidt  number  is 
assumed  to  be  a  constant  (  0.7  <  Scx<0.9  );  this  circumstance  allows  us  to  evaluate  the  turbulent  mass  transport  rate  by 
means  of  the  eddy  viscosity.  This  last,  in  turn,  is  a  local  flow  property  and  must  be  modelled  carefully  by  employing 
one  of  the  literature-known  models.  Following  the  Prandtl  mixing-length  theory  the  eddy  viscosity  was  assumed  to  be; 

Vt=IJt  (1.20) 
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where  l„  is  the  mixing  length  and  Vt  the  turbulent  velocity.  Typically  for  flows  near  a  wall  Vj  is  assumed  to  be 
proportional  to  the  product  of  the  mixing  length  and  the  magnitude  of  the  mean  velocity  gradient: 


Vr=l, 


du 

dy 


(1.21) 


The  formulation  of  both  l„  and  Vj  is  made  assuming  the  boundary  layer  to  be  a  composite  layer  consisting  of  inner  and 
outer  regions  merging  in  a  overlap  layer  which  represents  the  turbulent  boundary  layer  core.  The  extension  of  this  last 
region  increases  with  Re  along  the  boundary  layer  thickness. 

To  describe  the  velocity  and  gradients  profiles  through  the  boundary  layer  is  useful  to  define  viscous  scales  that  are  the 
appropriate  velocity  scales  and  length  scales  in  the  near-wall  region.  These  are  the  friction  velocity 


* 

u 


and  the  viscous  length  scale: 


Thus  the  distance  from  the  wall  measured  in  viscous  lengths  (or  wall  units)  is  defined  as: 


* 


(1.22) 


(1.23) 


(1.24) 


while  the  non  dimensional  local  mean  velocity  is: 


+  U 

U  =—  (1.25) 

u 

Inner  and  outer  regions,  together  with  their  sublayers,  are  defined  on  the  basis  of  y'^:  in  the  viscous  wall  region  < 
50)  shear  stresses  are  strongly  dependents  on  the  molecular  viscosity;  in  the  outer  layer  >  50)  the  viscosity  effect  is 
negligible. 

Within  the  viscous  wall  region,  the  Reynolds  shear  stress  is  negligible  compared  with  the  viscous  stress  in  the  viscous 
sublayer  (y"^  <  5):  they  are  of  the  same  order  in  the  buffer  layer  (5  <  y"^  <  50).  Furthermore  increasing  Re  the  fraction  of 
the  boundary  layer  occupied  by  the  viscous  wall  region  decreases:  in  fact,  y"^  being  similar  to  a  local  Reynolds  number, 
we  can  write: 
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1  =- 
^  lv=<5  S 


u*S 


=  Re 


* 


V 


(1.26) 


showing  that  S^,  j S  varies  as  Re*'*. 

Using  the  viscous  scales  together  with  equations  (IV. 24)  and  (IV. 25)  the  mixing  length  definitions  and  velocity 
gradients  obtained  by  the  universal  law  of  the  wall  can  be  summarized  in  each  sublayer  as: 


Inner  layer  -  viscous  (laminar)  sub-layer  (  5):  in  this  region  velocity  and  transversal  velocity  gradient  are 

given  by: 


and 


du 

dy 


*2 


(1.27) 


while  the  mixing  length  is  defined  by  the  van  Driest  equation: 


L  =  ky 


f 

l-exp 

V 


(1.28) 


where  the  suggested  constant  values  for  k  and  A'^  are  respectively  0.4  and  26. 


Inner  layer  -buffer  sub-layer  (  5  <3;^<  30)  :  in  this  region  velocity  and  transversal  velocity  gradient  are  given 
by: 


m'*' =  -3.05-1-5.0  In  and 


du 

dy 


0.2  j 


(1.29) 


while  the  van  Driest  mixing  length  equation  is  still  valid. 


Overlap  region  ( 30  <  ;  y  <  0.2  S ) :  here  velocity  and  transversal  velocity  gradient  are  given  by: 

du 


u^=  5.5 +  2.5  Iny" 


and 


dy 


0.4  y 


(1.30) 


and  the  mixing  length  is: 


L  =  ky 


(1.31) 


where  k  =  0.41  is  the  Kendall  constant. 
Outer  layer  -  wake  region  (y  >  0.2  ^): 
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=  5.5 +  —  In  + 
k 


n 

—  w 
k 


ySj 


(1.32) 


where  k  is  the  Kendall  constant  and  the  term  in  squared  brackets  is  a  function  called  the  law  of  the  wake;  in  particular 
for  zero  pressure  gradient  flow  the  wake  strength  parameter  11  =  0.48  and  the  wake  function  w  is: 


fy) 

1  •  2 

f  >) 

w 

=  2  sm 

[2S) 

(1.33) 


Thus  the  velocity  gradient  is: 


du 

u  .  2;r  n  .  1 

( TT  y'^ 

f ^  y") 

= - \-u - sm 

— 

cos 

— 

dy'_ 

k  y  S  k  ' 

^2S) 

while  the  mixing  length  is  assumed  to  be  a  constant  equal  to: 


I  =  0.09  S 

m 


(1.34) 


(1.35) 


Using  equations  (1.29)  -  (1.35)  together  with  (1.19),  (1.21)  velocity  profiles,  mixing  lengths  and  turbulent  mass 
diffusivities  were  calculated  along  the  boundary  layer  thickness  at  two  locations  (x  =  25  cm  and  x  =  50  cm);  the 
velocity  range  available  in  each  region  of  the  boundary  layer  at  x  =  25  is  summarized  in  table  4  while  the  turbulent  mass 
diffusivities  profiles  are  plotted  in  figure  6. 


Vert.  Distance  [mm] 

Velocity  [m/s] 

V 

+ 

y  <  0.025 

0  <  u  <  350 

5  <  <  30 

0.025  <  y  <0.135 

350  <  u  <  750 

V<  0.2  6 

0.135  <  y  <  1.8 

750  <  u  <  1100 

V  >  0.2  6 

1.8  <  y  <  9 

1100  <  u  <  Ve 

Table  4 
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O.E+OO  1.E-02  2.E-02  3.E-02  4.E-02  5.E-02  6.E-02  7.E-02 

Dt,  Turbulent  diffusivity  [m^/s] 


Fig  6:  Turbulent  mass  diffusivity  profiles  at  x  =  25  and  x  =  50 

As  shown  by  the  plot  the  turbulent  mass  diffusivity  is  more  than  one  order  of  magnitude  larger  than  the  hydrogen 
molecular  diffusion  coefficient;  maximum  values  are  localized  in  the  core  of  the  overlap  region  where  the  molecular 
viscosity  is  negligible.  Furthermore  Dj  increases  with  the  longitudinal  distance  allowing  diffusion  time  lower  than 
convective  time  along  all  the  grain. 


1.3  Liquid  Lithium  Droplets 

While  the  gaseous  hydrogen  diffuses  and  burns,  droplets  of  liquid  lithium  produced  by  LiH  thermal  decomposition  are 
dragged  away  along  the  combustion  chamber.  The  analysis  of  their  behaviour  is  important  to  understand  the  impact  of  a 
liquid  phase  on  the  engine  final  performance. 

In  fact  due  to  its  high  boiling  temperature  and  latent  heat  of  vaporization  (see  [2]),  liquid  lithium  acts  as  a  heat  sink: 
depending  on  thermal  fluxes,  droplet  diameter  and  time  of  residence  into  the  combustion  chamber,  it  can  exit  as  a 
unburnt  liquid  from  the  nozzle  taking  away  the  heat  trapped  along  the  path;  otherwise  it  can  vaporize  (partially  or 
completely)  and  burn  so  releasing  heat  useful  to  increase  the  flow  temperature. 

Moreover  from  the  point  of  view  of  the  mathematical  model  development,  the  presence  of  a  liquid  phase  evolving  into  a 
gaseous  stream  poses  nontrivial  problems  that  must  be  simplified  by  modelling  the  droplet  behaviour. 

A  brief  analysis  of  droplets  evolution  based  on  simplified  models  is  presented  in  the  following  sections. 

1.3.1  Droplets  Diameter 

Assuming  known  the  temperature  path  and  the  thermophysical  properties  of  the  heterogeneous  mixture,  the  driving 
factor  in  the  droplet  lifetime  analysis  is  the  droplet  diameter.  Large  drops  could  vaporize  and  burn  with  characteristic 
time  longer  than  the  particle  residence  time  in  the  combustion  chamber. 

To  asses  the  approximate  dimension  of  the  droplets  for  the  case  under  examination  we  move  from  some  simple 
considerations: 
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Depending  on  the  solid  grain  solution  design  we  can  decide  a  priori  the  most  effective  granulometry.  As  an 
example,  if  LiH  is  dispersed  into  HTPB  or  other  polymeric  binders,  we  can  use  (as  in  the  case  of  AP/HTPB  solid 
propellant)  LiH  particles  having  diameters  of  the  order  of  few  micrometers.  These  particles  once  liquefied, 
produce  droplets  with  diameters  of  the  same  order  of  magnitude. 

When  designed  without  a  binder,  the  whole  grain  surface  is  composed  of  solid  LiH:  as  shown  in  the  preceding 
sections  thermal  fluxes  liquefy  the  grain  surface  producing  a  thin  film  of  liquid  LiH  which  height  is  of  the  order 
of  few  micrometers.  When  dragged  away,  this  film  can  produce  droplets  with  maximum  diameters  of  the  same 
order  of  magnitude. 

As  shown  in  the  experiments  presented  (see  [2])  the  gaseous  hydrogen  formed  by  thermal  decomposition  of  LiH 
bubbles  violently  into  the  droplet,  favouring  its  fragmentation  into  smaller  particles. 


Moreover  these  droplets  when  dragged  away  from  the  surface  are  thrown  into  a  high  speed  turbulent  stream;  the  flow 
dynamic  pressure  acting  on  the  droplet  surface  contributes  to  its  fragmentation  into  smaller  droplets. 

The  problem  of  aerodynamic  break-up  of  liquid  drops  has  a  large  literature.  Most  works  analyze  the  break-up  dynamics 
as  a  function  of  the  ratio  between  the  aerodynamic  force  acting  on  the  particle  and  the  restoring  force  due  to  surface 
tension;  this  non  dimensional  ratio  is  defined  as  the  Weber’s  number  : 


We, 


(1.36) 


where  D,  v,  o  and  p  are  respectively  the  droplet  diameter,  velocity  (relative  to  the  stream),  surface  tension  and  the  air 
density.  As  stated  by  Pilch  &  Erdman  (1987)  the  break-up  dynamic  by  means  of  shadowgraph  techniques, increasing 
the  flow  speed  the  break-up  dynamic  range  from  the  simple  separation  into  smaller  droplets  (We  ==12)  up  to  a 
catastrophic  nebulization  into  tens  of  very  small  droplets  (We  >  350);  conversely,  for  assigned  liquid  properties  and 
flow  speed  there  is  only  one  range  of  drop  diameters  compatible  with  a  specific  break-up  dynamics  regime.  It  is 
possible  thus  to  define  a  Weber  number  (called  Critical  Weber’s  Number  We*^ )  under  which  the  drop  break-up  cannot 
be  observed  within  times  significant  in  practical  applications.  An  experimental  correlation  between  the  critical  Weber 
and  the  Ohnesorge  number:’^ 


Oh= 


1^1 


VPi 


(1.37) 


where  |j.i,  pi  and  D  are  respectively  the  droplet  viscosity,  density  and  diameter  was  proposed  by  Brodkey  (Figure  1) 


17 


Grant  Number  FA8655-10-1-3091 


Fig.  7:  Brodkey  correlation 


As  shown  in  the  picture,  the  critical  Weber  is  about  12  for  Ohnesorge  numbers  <0.1,  increasing  dramatically  for  Oh  > 
0.1.  In  other  words,  if  Oh  >0.1  (e.g.  when  the  fluid  viscosity  increases)  the  drops  break-up  becomes  quite  impossible. 
The  Brodkey  correlation  can  be  expressed  analytically  as: 


We;;  =  12  (1-1-1.077  Oh'®) 


(1.38) 


To  determine  the  maximum  diameter  compatible  with  the  flow  condition  for  the  case  under  examination,  we  have 
calculated  using  (1.37)  and  (1.38)  the  critical  liquid  lithium  Ohnesorge  number  as  a  function  of  the  drop  diameter, 
assuming  its  viscosity,  density  and  surface  tension  at  the  LiH  decomposition  temperature  (about  1100  K);  then  the  flow 
Weber  number  was  calculated  as  a  function  of  the  drop  diameter  for  air  speeds  ranging  from  500  up  to  1500  m/s 
assuming  its  density  valued  at  the  air  stream  temperature  at  the  combustion  chamber  entrance  (T  =1400  K).  Results  are 
plotted  in  figure  8. 


Fig  8:  Liquid  Li  drop  diameter  compatible  with  the  SCRJ  flow  conditions  assumed 
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As  shown  in  the  picture  for  high  speed  flow  the  Weber  number  of  droplets  having  diameters  >  50  pm  is  always  greater 
than  the  critical  Weber  number;  in  particular,  for  the  flow  conditions  expected  in  the  core  of  the  boundary  layer,  the 
maximum  allowed  droplet  diameter  is  of  order  of  10  pm. 


I. 3.2  Droplet  evaporation  and  burning 

Once  the  droplets  maximum  diameter  is  known  we  can  estimate  its  gasification  rates  and  lifetime,  a  key  issue  in 
evaluating  the  heat  subtracted  or  added  by  lithium  to  the  flow  while  moving  along  the  combustion  chamber. 
Incomplete  evaporation  could  in  fact  signify  potential  heat  trapped  into  the  drop  and  transported  out  of  the 
combustor,  thus  ineffective  for  propulsion.  Conversely,  fast  evaporation  yields  gaseous  lithium  that  can  react  while 
vaporizing,  or  diffuse  away  from  the  drop  (e.g.;  when  the  oxidizer(s)  concentration  is  too  small  to  allow  any 
reaction)  and  burn.  This  second  case  allows  us  to  justify  the  idea  of  lithium  hydride  as  a  powerful  dual-fuel  system: 
in  fact  we  can  roughly  consider  LiH  as  a  special  hydrogen  storage  system  in  which  the  storage  medium  (lithium)  is, 
in  turn,  a  highly  energetic  fuel. 

With  appropriate  and  simple  assumptions,  the  evaporation  and  burning  of  spherical  liquid  droplets  is  relatively 
simple  to  analyze;  for  both  cases  we  have  used  models  based  on  closed-form  analytic  solutions  of  the  simplified 
governing  conservation  equations  and  that  are  well  know  in  the  literature.^*’  These  solutions  have  allowed 
us  to  explore  the  influence  of  droplet  size  and  ambient  conditions  on  droplet  evaporation  and  burning  times  (see 
Appendix  Al)  and  to  definitively  conclude  that  almost  all  liquid  drops  vaporize  and  burn  before  exhausting  from  the 
combustion  chamber,  and  contribute  to  the  heat  release  in  the  airstream. 

II  Numerical  simulation  -  Introduction 

As  stated  in  [2],  the  attention  of  almost  all  supersonic  combustion  researchers  being  focused  on  liquid  fuels,  numerical 
and  experimental  data  related  to  solid  fuel  scramjet  study  are  quite  scarce  or  nonexistent  in  the  open  literature.  For  this 
reason,  to  simulate  the  LiH  fuelled  SFSCRJ  we  have  investigated  the  geometrical  configurations  and  the  “flight” 
conditions  presented  in  [1].^*'^^ 

In  this  section  we  will  focus  on  issues  associated  with  the  boundary  conditions,  with  particular  attention  to  the  injection 
of  mass  from  the  grain  surface  (blowing  effect). 

Then  the  numerical  method  will  be  described,  stressing  the  solution  adopted  (the  “sources”  method)  to  account  for  the 
blowing  and  its  effect  on  the  boundary  layer. 

II. l  Grain  regression  model 

From  the  physical  viewpoint  a  solid  fuelled  scramjet  presents  many  similarities  with  hybrid  rockets;  the  oxidizer 
diffuses  into  the  boundary  layer  and  reacts  with  the  gaseous  fuel  pyrolyzed  and  released  by  the  solid  grain  surface.  It  is 
thus  reasonable  to  assume  (accounting  for  the  specific  differences)  one  among  the  well  known  models  present  in  the 
open  literature  to  investigate  the  surface  regression. 
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The  most  plausible  (and  used)  model  of  hybrid  combustion  was  developed  by  Marxman  and  Gilbert  in  1963.^®’  It  is 
based  on  the  concept  of  diffusive  flame;  in  a  hybrid  system  the  combustion  is  localized  in  a  narrow  zone  of  the 
boundary  layer  where  the  air  oxygen  diffusing  from  the  main  stream  meets  the  fuel  released  by  the  surface.  As  stated  by 
Carmicino  this  feature  allows  to  study  of  flame  zone  as  a  discontinuity  of  the  temperature  gradient  and  of  the  mixture 
composition.®** 

In  this  model  the  grain  regression  rate  is  governed  by  the  local  thermal  flux  at  the  surface;  in  fact  we  can  write: 


where  Pg  is  the  density  and  Vg  the  velocity  of  the  gaseous  fuel  leaving  the  surface;  (J/m^s)  is  the  thermal  flux  at  the 

wall  while  (J/kg)  is,  in  general,  the  total  heat  required  to  vaporize  the  solid  fuel.  In  the  case  of  a  LiH  grain,  the 
surface  temperature  being  that  of  the  phase  change  (LiH  solid  to  liquid;  950  K),  accounts  for  the  heat  necessary  to 
vaporize  and  decompose  LiH:  as  discussed  in  [2],  //„=!!  MJ/kg. 

The  equation  describing  the  regression  rate  r  is  obtainable  using  equation  IT  1 : 

=  ("2) 

that  is: 


r  =  (11.3) 

Ps  ti. 

Thus  the  key  point  of  the  model  is  the  evaluation  of  thermal  flux  at  the  surface;  neglecting  the  conduction  heat 
exchange  inside  the  grain,  total  thermal  flux  is  the  sum  of  radiative  and  convective  heat  fluxes.  This  last,  in  turn,  is 
affected  by  the  injection  of  gaseous  fuel  (at  950  K)  from  the  surface. 


II.2  Blowing 

To  assess  the  impact  of  the  fuel  injection  (also  known  as  blowing)  from  the  surface  on  the  convective  heat  flux,  we 
must  account  for  its  influence  on  the  local  velocity  profile  in  the  boundary  layer.  In  fact  the  heat  transfer  coefficient  is 
determined  using  the  Stanton  number,  which,  in  turn,  is  computed  from  the  skin  friction  coefficient  calculated  by  using 
the  wall  functions. 

The  first  theoretical  analysis  of  blowing  was  made  in  1942  by  H.  Schlichting;  successively  experimental  as  well  as 
theoretical  investigation  have  been  performed  by  J.  C.  Rotta.  Figures  9-10  show  some  of  Rotta’s  results. 

Injecting  gas  from  the  surface,  the  boundary  layer  thickness  grows.  Figure  V.4  shows  the  variation  of  momentum 
thickness  82(x)  along  a  porous  flat  plate  with  homogeneous  blowing  (and  suction)  at  various  values  of  the  blowing 
velocity  (v,^)  at  the  wall.  In  this  case  the  external  velocity  was  Ue  =  20  to  30  m/s  and  the  normal  wall  velocity  ranged 
from  v„  =  -0.1  m/s  (suction)  up  to  0.13  m/s  (blowing). 
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As  a  consequence  of  this  variation,  the  velocity  profiles  in  the  boundary  layer  change;  to  calculate  the  new  profiles, 
Rotta  performed  a  theoretical  formulation,  based  on  the  Prandtl’s  mixing  length  assumption,  that  is  summarized  in  the 
following  equation:^* 


u  I  1  V  r  1  I 

—  =  -  Inj]  +C+ - ^  -  \nr/  +  C  (II.4) 

V,  ^  4  \x 


where  JJ  =  y  W^jv  is  the  nondimensional  distance  from  the  wall,  v  the  viscosity,  and  V,  =  p)'^  is  the  friction 

velocity. 

Curves  calculated  by  Rotta  using  this  equation  are  plotted  in  figure  10  , while  figure  1 1  shows  a  comparison  between 
theoretical  and  experimental  data. 
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Fig.  10  :  Velocity  distribution  in  the  B.L.  according  to  eqn.  II.IO 


L0.2;28.2  100^ 


Fig  11 :  Comparison  between  theoretical  (lines)  and  experimental  (circles)  data 


The  injection  of  a  gas  through  a  porous  wall  into  a  compressible  flow,  up  to  Mach  3.6,  was  investigated  by  L.  C.  Squire; 
its  calculations  showed  that  Rotta’s  assumption  were  still  valid  leading  to  satisfactory  results. 

So,  what  happens  really  in  a  scramjet  combustor  chamber?  High  temperature  (1450  K)  and  high  flow  velocity  (1400 
m/s)  indicate  high  friction  velocities  along  the  grain  (up  to  90  m/s).  At  the  same  time,  for  regression  rates  of  order  1 
mm/s  (or  less),  from  equation  II.4  we  can  see  that  the  expected  blowing  velocities  are  of  order  of  few  m/s:  thus  the 
V^/v*  ratio  is  very  small,  of  order  of  0.01.  Now,  observing  figure  1 1 ,  we  can  see  that,  for  T]  <  20,  the  w/v*  ratio  is  < 

0.1.  This  means  that  the  blowing  impact  on  velocity  profiles  and  heat  fluxes  at  the  wall  is  (in  this  case!)  negligible.  This 
assumption  however  does  not  account  for  the  cooling  effect  due  to  the  gas  injection  temperature;  furthermore  there  are 
zones  in  which  the  boundary  layer  behaviour  does  not  agree  with  the  theoretical  model  developed  for  a  flat  plate  (e.g:  in 
the  flameholder  recirculation  zone  or  in  regions  with  shock-boundary  layer  interactions).  In  these  regions,  the  blowing 
effect  will  be  investigated  analyzing  simulation  results. 


II.3  Numerical  approach 

Our  simulations  will  be  carried  on  using  FANS  and  an  axisymmetric  domain.  Numerical  computations  are  carried  out 
using  the  Fluent'’''^  6.3  computational  fluid  dynamics  package. 
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II.3.1  Geometries 

The  geometries  considered  in  this  analysis  (hereinafter  referred  respectively  as  cylindrical  and  dump  combustor)  are 
shown  in  figure  12  and  13. 

Figure  12  shows  the  longitudinal  half-section  of  the  cylindrical  combustor. 


Fig  12:  Cylindrical  combustor 

As  in  the  cases  investigated  by  Jakymowycz,  Yang  and  Kuo,  and  by  Ben-Arosh  and  Gany,  the  LiH  grain  is  located 
along  the  upper  surface  of  the  2-D  half-section,  with  an  inlet  height  of  5  cm;  the  LiH  grain,  60  cm  long,  is  located  after 
a  short  adiabatic  wall  (5  cm,  the  inlet  prosecution);  then  another  adiabatic  wall  (or  inert  plate),  spanning  40  cm  along 
the  chamber,  is  the  final  part  of  the  combustor. 

Figure  13  shows  the  longitudinal  half-section  of  the  dumped  combustor  (configuration  Dumpl). 


5  cm  60  cm  40  cm 


Fig  13:  Dump  combustor  (coufiguratiou  Dl) 

The  dump  combustor  differs  from  the  cylindrical  one  for  the  presence  of  a  sudden  expansion,  with  a  dump  height  of  1 
cm. 

We  have  also  investigated  a  third  dump  configuration  (configuration  DumpOS),  that  differs  from  Dumpl  only  in  the 
dump  height  (0.5  cm). 

The  numerical  technique  is  a  finite  volume  approach  with  quadrilateral  control  volumes  and  structured  mesh.  For  all  the 
calculations  we  use  a  Cartesian  grid.  In  particular  we  have  used  three  different  meshes: 

The  cylinder  configuration  (fig  12)  has  a  grid  composed  by  155000  cells:  in  the  longitudinal  direction  over  the  grain  the 
grid  spacing  is  regular  with  a  uniform  0.5  mm  mesh.  The  mesh  increases  downstream  the  grain,  with  the  largest  grid 
measuring  2  mm  in  length.  In  the  vertical  direction,  the  grid  has  an  exponential  progression  with  an  exponential  ratio  of 
0.17  and  the  smallest  grid  (near  the  surface)  measuring  50)J,m  in  height. 
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The  Dump  1  configuration  of  the  dump  combustor  (fig  13)  has  a  grid  composed  by  216900  cells:  in  the  longitudinal 
direction  over  the  grain  the  grid  spacing  is  regular  with  a  uniform  0.25  mm  mesh.  The  mesh  increases  downstream  the 
grain,  with  the  largest  grid  measuring  2  mm  in  length.  In  the  vertical  direction,  the  grid  has  an  exponential  progression 
with  an  exponential  ratio  of  0.1 1. 

Finally  the  Dump05  configuration  of  the  dump  combustor  has  a  grid  composed  by  162900  cells:  in  the  longitudinal 
direction  over  the  grain  the  grid  spacing  is  regular  with  a  fixed  mesh  of  0.25  mm.  The  mesh  increases  downstream  the 
grain,  with  the  largest  grid  measuring  2  mm  in  length.  In  the  vertical  direction,  the  grid  has  an  exponential  progression 
with  an  exponential  ratio  of  0.1 1. 

Even  though  a  flow  with  injection  of  fluid  normal  to  a  wall  can  lead  to  intrinsic  instabilities,**  we  have  assumed  a 
steady  state  approach  using  a  coupled,  implicit,  second  order  upwind  formulation  scheme.  The  turbulent  equations  are 
modelled  using  a  standard  k-e  model.  ’  ’ 

The  method  adopted  to  implement  the  boundary  condition  at  the  grain  surface  requires  some  detailed  considerations  and 
will  be  presented  in  the  following  section. 

For  the  gas  phase  combustion  we  have  adopted  a  six  species  simplified  mechanism  involving  the  following  reactions: 
LiH  decomposition:"*^’ 


LiH  ^  Li  +  1/2  Hz 


(11.5) 


One-step  reaction  for  hydrogen  combustion:^^’ 


Hz  H- 1/2  Oz  ^  HzO 


(11.6) 


and  one-step  reaction  for  lithium  oxidation:^*’ 


2Li  H- 1/2  Oz  ^  LizO  (II.7) 

This  mechanism  does  not  account  for  backward  reactions  since  their  characteristic  times  (as  stated  in  [2])  are  higher 
than  flow  convective  time.  The  coupling  between  chemical  kinetics  and  turbulence  is  modelled  using  the  Eddy 
Dissipation  Concept  (EDC)  model.  A  description  of  the  Fluent'’''^  solver  among  with  k  -e  and  EDC  models  is  available 
in  Fluent  User’s  Guide  manual.*^ 

Thermodynamic  data  (formation  enthalpy  and  entropy  at  298.15  K,  constant  pressure  specific  heat)  data  for  each 
chemical  species  have  been  calculated  using  the  CEA2  database;  in  particular,  we  have  calculated  the  constant  pressure 
specific  heats  from  300  K  up  to  6000  K  and  successively  we  have  implemented  them  in  Fluent  by  means  of 
polynomials. 

In  the  same  range  of  temperatures  (300K  -6000K)  we  have  calculated  all  the  molecular  viscosities  p  using  the  following 
equation  (from  kinetic  theory): 

2  1  /mkT 
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where  k=  1.38  x  10'^^  [J/K]  is  the  Boltzmann’s  constant,  m  the  mass  of  a  molecule  and  o  =  Ttd^ ,  with  d  molecular 
diameter.  Curves  obtained  have  been  thus  interpolated  using  third  order  polynomial  and  then  implemented  in  Fluent. 
Finally  we  have  calculated  for  all  species  the  binary  molecular  diffusion  coefficients  from  300  K  to  6000K  using  the 
model  presented  in  1.3.2. 1  (equations  1.16  -  1.18).  These  values  have  been  also  interpolated  by  means  of  third  order 
polynomial  and  implemented  in  Fluent  as  well. 


II.3.2  Blowing  model 

In  order  to  model  injection  of  fluid  from  the  wall,  we  have  to  account  for  the  characteristic  length  scales  involved  in  the 
process.  If  the  characteristic  length  between  two  injection  points  is  much  smaller  than  the  length  scales  of  the  turbulent 
channel  flow,  the  details  of  injection  at  this  scale  are  not  important  and  the  injection  can  be  modelled  by  a  continuous 
model.  Furthermore,  the  model  must  account  for  the  viscous  drag,  together  with  the  mass,  momentum  and  enthalpy 
generation  from  the  fluid  injected.  In  our  study,  we  have  modelled  directly  the  physical  phenomenon,  instead  of  using 
additional  terms  in  the  wall  functions. 

The  blowing  model  adopted  ,  the  so-called  ‘sources  model’,  consists  in  applying  mass,  momentum  and  energy  sources 
at  the  first  cell  above  an  impermeable  adiabatic  wall 

We  have  implemented  the  sources  terms  writing  in  C  language  an  UDF  (User  Defined  Function).  A  user-defined 
function,  or  UDF,  is  a  function  that  we  program  and  that  can  be  dynamically  loaded  with  the  solver  to  change  or 
enhance  the  standard  features  of  the  code.  Our  UDF  operates  on  the  source  term  of  conservation  equations;  that  is, 
writing  these  equations  in  a  generalized  form  we  have: 

div[/7y(^-r^  grad(/>]  =  5^  (II.8) 

where  the  first  and  second  left-hand  side  term  represent  convection  and  diffusion,  respectively,  and  the  right-hand  side 
term  is  a  source  term.  In  the  equation  (/)  stands  for  the  conserved  property  and  the  appropriate  diffusivity  coefficient. 
To  create  an  appropriate  model  of  blowing,  we  have  defined,  among  with  a  mass  source  (for  every  cell  over  the  grain 
surface),  also  a  momentum  source,  a  species  source  and  an  energy  source.  In  fact,  defining  only  a  mass  source,  that 
mass  would  enter  the  domain  with  no  momentum  or  thermal  energy.  The  mass  will  therefore  have  to  be  accelerated  and 
heated  at  the  expense  of  the  main  airstream  into  which  is  injected,  and  consequently  there  might  be  a  drop  in  velocity 
and/or  temperature.  This  drop  may  or  may  not  be  significant,  depending  on  the  size  of  the  source;  in  both  cases  it  could 
modify  significantly  the  solution.  Finally,  the  energy  conservation  being  written  in  terms  of  total  energy,  the  energy 
source  term  must  account  also  for  the  kinetic  energy  of  the  gas  leaving  the  surface. 

The  source  term  adopted  for  the  variables  are: 

Mass  source  — 

n  2  A 

Y  -  Momentum  source  — 
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Species  (LiH)  source 


Energy  source  : 


S 


energy 


Pg'^gKso 


where  A  and  V  are,  respectively,  the  area  of  the  cell  surface  adjacent  to  the  solid  grain  and  the  cell  volume. 

These  sources  are  directly  connected  with  equation  II. 1  and  II. 3.  Using  the  equation  II. 1  we  define  the  gas  mass  flow 
rate  as  a  function  of  the  heat  flux  arriving  on  the  surface.  Analogously,  using  equation  II. 3,  we  can  calculate  the  local 
regression  rate  as  a  function  of  the  surface  heat  flux.  The  key  point  is  thus  how  to  calculate  the  heat  fluxes  on  the  basis 
of  velocity,  density,  pressure  and  temperature  field  calculated  by  the  solver  in  the  gas  above  the  surface. 

To  accomplish  this  task,  we  have  implemented  in  the  UDFs  the  same  model  (based  on  the  standard  wall  function)  used 
by  the  solver  to  calculate  the  total  surface  heat  flux  and  summarized  in  the  following. 

The  wall  functions  in  the  code  are  based  on  the  proposal  of  Launder  and  Spalding:*^  the  law-of -the- wall  for  the  mean 
velocity  yields 


where 


with; 


U-='-  ln(£/) 


1 1 

u  =— — ^ 

^Jp 

P  cf  kf  yp 


y  =■ 


p 


(11.9) 


(II.  10) 


(11.11) 


k 

von  Karman  constant  (=  0.4187) 

E 

empirical  constant  (=  9.793) 

Up 

mean  velocity  of  the  fluid  at  point  P 

kp 

turbulence  kinetic  energy  at  point  P 

yp 

distance  from  point  P  to  the  wall 

p 

dynamic  viscosity  of  the  fluid 

Cm 

k-e  model  constant  (=  ) 
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The  logarithmic  law  for  mean  velocity  is  known  to  be  valid  for  30  <  y*  <300.  In  the  code  ,  the  log-law  is  employed 
when  y*  >  1 1.225.  When  the  mesh  is  such  that  y*  <  1 1.225  at  the  wall-adjacent  cells,  Fluent  applies  the  linear  (laminar) 
stress-strain  relationship  that  can  be  written  as 


*  * 

U  =y  (11.12) 

We  write  a  similar  logarithmic  law  for  mean  temperature  by  using  the  Reynolds  analogy  between  momentum  and 
energy  transport;  as  in  the  law-of-the-wall  for  the  mean  velocity,  the  wall  functions  for  temperature  employed  by  Fluent 
comprises  the  following  two  different  laws: 


•  linear  law  for  the  thermal  conduction  sublayer,  where  molecular  conduction  is  important 

•  logarithmic  law  for  the  turbulent  region  where  the  effect  of  turbulence  dominate. 

In  highly  compressible  flows,  the  temperature  distribution  in  the  near-wall  region  can  be  significantly  different  from 
that  of  subsonic  flows,  due  to  the  heating  by  viscous  dissipation.  In  Fluent,  the  temperature  wall  functions  include  the 
contribution  also  from  the  viscous  heating  and  have  the  following  composite  form: 


T  = 


(T„-T,)pcCl‘kf 


(11.13) 


1 

r=Pr/+-/?Pr  ^  ^  Ul 

2  q 


/  *  *  \ 

(t  <Tr)  (n-14) 


T  =Pr 


iln(£/)  +  P 


+  i/7jy^{Pr,C/;+(Pr-Pr,)(/*}  (y  > /r) 


(11.15) 


where  P  is  computed  by  using  the  formula  given  by  Jayatilleke: 


P 


9.24 


1  [l  +  0.28e“'’“"’’'/'''''] 


(11.16) 


Prt 

turbulent  Prandtl  number  (0.85  at  the  wall) 

Pr 

Molecular  Prandtl  number  (|a,  Cp/  kf) 

Cp 

specific  heat  of  fluid 

27 


Grant  Number  FA8655-10-1-3091 


Uc 

mean  velocity  magnitude  at  y*  = 

kp 

turbulent  kinetic  energy  at  point  P 

yp 

distance  from  point  P  to  the  wall 

M 

dynamic  viscosity  of  the  fluid 

Cp 

k-e  model  constant  (=  ) 

Tp 

temperature  at  the  cell  adjacent  to  wall 

T„ 

temperature  at  the  wall 

P 

density  of  fluid 

*  * 

The  non-dimensional  thermal  sublayer  thickness,  ,  in  equations  11.14,15  is  computed  as  the  y  value  at  which  the 

linear  law  and  the  logarithmic  law  intersect,  once  the  molecular  Prandtl  number  of  the  fluid  being  modelled  is  known  or 
assumed. 

The  procedure  used  by  Fluent  to  apply  the  law-of-the-wall  for  temperature  is  as  follows.  Once  the  physical  properties  of 
the  fluid  are  specified,  its  molecular  Prandtl  number  is  calculated.  Then,  given  the  molecular  Prandtl  number,  the 
thermal  sublayer  thickness  is  computed  from  the  intersection  of  the  linear  and  logarithmic  profiles,  and  stored.  During 
the  iteration,  depending  on  the  y*  value  at  the  near-wall  cell,  either  the  linear  or  the  logarithmic  profile  in  Equation 
11.14  or  11.15  is  applied  to  compute  the  heat  flux  (depending  on  the  type  of  the  boundary  conditions  on  the  wall). 


II.3.3  Test  Case 

To  simulate  the  LiH  fuelled  SFSCRJ  and  to  compare  our  results,  we  have  investigated  geometrical  configurations  and 
“flight”  conditions  similar  to  those  presented  in  [1]  in  particular  we  have  selected  as  inlet  conditions  the  following:: 


Inlet 

Static  Temperature 

Static  Pressure 

Total  Pressure 

Flow  velocity  [m/s] 

Mach 

[K] 

[atm] 

[atm] 

2 

1550 

1 

7.8 

1490 

Table  6:  Inlet  conditions 


These  conditions  have  been  applied  to  each  geometrical  configuration  presented  in  this  section. 

III.  Results 

Figures  14  and  15  show  the  schematic  geometry  and  the  reference  system  used  (here  reported  for  the  cylindrical  and  for 
the  Dumpl  configuration). 
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Fig.  14:  Cylindrical  configuration:  reference  system 
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Fig.  15:  Dumpl  configuration:  reference  system 

The  fuel  grain  (red  line)  spans  60  cm,  from  x  =  5  to  x=  65  cm;  the  dump  ratio  Dr,  is  defined  as: 


D 


R 

_ 

R 

1 


where  Rc  is  the  combustor  internal  radius  and  Ri  is  the  inlet  internal  radius. 

Furthermore  we  have  defined  a  control  line,  (half  radius  line)  located  at  Rc/2,  that  we  will  use  to  plot  pressure  and 
temperature  along  the  combustor. 

Finally,  for  all  the  configurations,  we  have  performed  comparisons  between  reacting  and  non  reacting  cases,  assuming 
the  non  reacting  case  to  represent  the  flow  condition  immediately  before  the  start  of  the  surface  decomposition. 


111.1  Cylindrical  configuration  (Dr  =1) 


Figure  16  shows  the  static  temperature  contours  of  both  non  reacting  and  reacting  case. 
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Fig.  16:  Contours  of  static  temperature  (reacting  and  non  reacting  case) 


As  we  can  see,  an  intense  flame  is  present  over  the  fuel  grain,  growing  along  the  combustor  and  up  to  its  end.  As  an 
effect  of  the  gas  phase  combustion,  the  heat  flux  (figure  17)  is  much  higher  (up  to  three  times)  than  that  due  only  to  the 
hot  airstrem  entering  the  combustor  in  the  non  reacting  case. 


Fig.  17:  Heat  flux  over  the  grain  surface 


As  a  consequence  of  this  significant  heat  flux,  the  grain  surface  decomposes  and  regresses  with  a  regression  rate  of  the 
same  order  of  that  calculated  by  Jarymowycz  and  Kuo;^*  figures  18  -  20  show  respectively  the  regression  rate,  the  LiH 
mass  flow  rate  produced  and  the  associated  blowing  velocity. 
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Fig.  18:  Grain  regression  rate 


Fig.  19:  LiH  local  mass  flow  rate 
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Fig.  20:  Blowing  velocity 


The  total  LiH  mass  flow  rate  injected  into  the  chamber  is  0.065  kg/s  while  the  mass  flow  rate  of  the  air  entering  the 
combustor  is  2.77  kg/s;  thus  the  scramjet  is  working  at  an  air  to  fuel  mass  ratio  (O/F)  of  42.2 

The  maximum  temperature  reached  along  the  combustor  is  2950  K,  located  downstream  of  the  entrance  and  toward  the 
end  of  the  grain  (figure  21). 


5  15  25  35  45  55  65  75  85  95 

X  -  coordinate,  [cm] 


Fig.  21:  Maximum  temperature  along  the  chamber 


The  net  temperature  increase  for  x  >  65  (the  final  part  of  the  grain)  is  connected  to  the  end  of  the  cooling  effect  of  the 
gas  (which  is  at  a  temperature  of  950  K)  injected  in  the  reacting  boundary  layer  from  the  decomposing  and  regressing 
surface.  Furthermore,  as  a  consequence  of  the  ceased  mass  injection,  the  thickness  of  the  reacting  zone  decreases  slowly 
(figure  16),  allowing  an  expansion  of  the  core  flow.  This  effect  is  also  visible  in  figures  22,  23  where  we  have  reported 
the  temperature  calculated  along  the  chamber  axis  and  along  the  half  radius  line,  compared  with  the  non  reacting  case 
(that  is:  when  the  blowing  is  not  present),. 
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Fig.  22:  Static  Temperature  along  the  chamber  axis 
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Fig.  23:  Static  Temperature  along  the  half  radius  line 


Downstream  of  this  first  expansion,  the  reacting  zone  thickness  increases  due  to  the  increasing  combustion  temperature, 
allowing  a  certain  amount  of  flow  recompression;  this  effect  is  shown  in  figures  24,  (where  the  non  reacting  and  the 
reacting  case  are  compared),  and  in  figures  26,  27  (where  the  pressure  along  the  axis  and  the  half  radius  line  is  plotted). 
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Fig.  24:  Contours  of  static  pressure  (reacting  and  non  reacting  case) 


Fig.  25:  Contours  of  static  pressure  (detail  of  the  combustor  half  section) 


Figures  26  and  27  show  also  the  effect  of  the  complex  oblique  shock  system  on  initial  part  of  the  grain  (figures  25).  As 
we  can  see  an  oblique  shock  is  located  exactly  at  the  grain  inlet;  the  consequent  compression,  increases  as  an  effect  of 
the  waves  reflecting  along  the  chamber,  causing  the  decreasing  trend  of  the  blowing  velocity. 
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Fig.  26:  Static  pressure  along  the  chamber  axis 
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Fig.  27:  Static  pressure  along  the  half  radius  line 


This  feature  is  also  visible  in  figure  28  where  the  Mach  number  contours  for  the  non  reacting  and  reacting  cases  are 
compared.  Notice  that  despite  the  complex  wave  system  the  flow  at  the  exit  section  is  largely  supersonic. 
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Fig.  28:  Contour  of  Mach  number  (reacting  and  non  reacting  case) 


Finally,  figure  29  shows  the  contours  of  mixing  time,  calculated  as: 


t  = 


V 


2 


Dth 


where  Dth 

is  the  turbulent  diffusivity  and  Tj  is  the  Taylor  microscale,  calculated  as: 


v  = 


whit  V  kinematic  viscosity,  k  turbulent  kinetic  energy  and  e  the  turbulent  dissipation  rate. 
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Fig.  29:  Contour  of  mixing  time  (reacting  and  non  reacting  case) 


As  shown,  due  to  the  highly  turbulent  flow,  in  the  boundary  layer,  the  mixing  time  is  very  short  (of  order  of  10'^  s); 
moreover,  in  presence  of  combustion,  thanks  to  the  growth  of  the  boundary  layer  and  to  the  consequent  complex  wave 
system,  the  mixing  time  decreases  also  in  the  core  flow. 
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III.2  Dump  combustor  (configuration  DumpOS,  Dr  =  1.1) 


Figure  30  shows  the  static  temperature  contours  of  both  non  reacting  and  reacting  cases. 
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Fig.  30;  Contours  of  static  temperature  (reacting  and  non  reacting  case) 


As  we  can  see,  also  in  this  case  an  intense  flame  is  present  over  the  fuel  grain,  growing  from  the  recirculation  zone  near 
the  sudden  area  increase  up  to  the  exit  section.  As  an  effect  of  the  gas  phase  combustion,  the  heat  flux  (figure  31)  is 
much  higher  than  that  due  only  to  the  hot  flow  entering  the  combustor  (that  is,  in  the  non  reacting  case). 
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Fig.  31:  Heat  flux  over  the  grain  surface 


As  a  consequence  of  the  heat  flux,  the  grain  surface  decomposes  and  regresses  with  a  regression  rate  of  the  same  order 
of  calculated  by  Ben-Harosh  and  Gany;^^  figures  32  -  33  show  respectively  the  regression  rate,  the  LiH  mass  flow  rate 
produced  and  the  consequent  blowing  velocity. 


Fig.  32:  Local  regression  rate 
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Fig.  33:  LiH  local  mass  flow  rate 


Fig.  34:  LiH  local  mass  flow  rate 

The  total  LiH  mass  flow  rate  injected  inside  the  chamber  is  0.075  kg/s  while  the  mass  flow  rate  of  the  air  entering  the 
combustor  is  2.77  kg/s;  thus  the  scramjet  is  working  at  a  air  to  fuel  mass  ratio  (O/F)  of  37.05. 

The  maximum  temperature  reached  along  the  combustor  is  2906  K,  located  downstream,  near  the  end  of  the  grain 
(figure  35). 
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Fig.  35:  Maximum  temperature  along  the  chamber 

Once  again,  the  net  temperature  increase  after  x  =  65  (the  final  part  of  the  grain)  is  connected  to  the  end  of  the  cooling 
effect  of  the  gas  (which  as  a  temperature  of  950  K)  injected  in  the  reacting  boundary  layer  from  the  surface. 
Furthermore,  in  the  first  part  of  the  grain,  from  about  x  =  10  to  x=  20  cm,  the  cooling  effect  due  to  the  expansion  wave 
located  in  the  upper  corner  of  the  dump  is  well  visible.  This  effect  is  clearer  examining  figures  36,  37  where  we  have 
reported  the  temperature  calculated  along  the  chamber  axis  and  along  the  half  radius  line,  compared  with  the  non 
reacting  case. 
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In  fact,  as  we  can  see  comparing  figures  30  and  36  -  38,  this  expansion  wave  lowers  the  core  flow  temperature  up  to 
1300  K;  the  following  reflecting  waves  (figure  38)  create  alternatively  compression  and  expansion  regions,  affecting 
greatly  the  reacting  zone.  This  feature  is  shown  more  clearly  in  figures  39,  40,  where  a  detail  of  the  first  part  of  the 
grain  is  stressed. 
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Fig.  38:  Contours  of  static  pressure  (reacting  and  non  reacting  case) 
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Fig.  39:  Contours  of  static  pressure  (detail  of  reacting  case) 


Fig.  40:  Contours  of  static  temperature  (detail  of  reacting  case) 


The  pressure  behaviour  along  the  axis  chamber  and  the  half  radius  line  is  shown  in  figures  41,  42. 
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Fig.  41:  Static  pressure  along  the  chamber  axis 
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Fig.  42:  Static  pressure  along  the  half  radius  line 


Finally,  figures  43  and  44  show  respectively  the  Mach  number  contour  and  the  mixing  time  contour  of  both  non 
reacting  and  reacting  cases. 
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Fig.  43:  Contour  of  Mach  number  (reacting  and  non  reacting  case) 


Observing  figure  43  we  can  see  the  effect  of  the  temperature  increase;  due  to  the  good  mixing  times  in  the  chamber 
(figure  44)  and  to  the  presence  of  a  recirculation  zone  at  the  entrance  inside  the  grain  we  obtain  higher  combustion 
temperatures  than  in  the  case  of  the  cylindrical  configuration,  while  maintaining  a  supersonic  flow  in  the  combustor  exit 
section. 
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Fig.  44:  Contour  of  mixing  time  (reacting  and  non  reacting  case) 
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III.3  Configuration  Dumpl  (  Dr  =  1.2)  and  comparisons 


Figure  45  shows  the  static  temperature  contours  of  both  non  reacting  and  reacting  case. 


Fig.  45:  Contours  of  static  temperature  (reacting  and  non  reacting  case) 


As  in  the  case  of  the  first  two  configuration  examined,  a  flame  is  present  over  the  fuel  grain,  growing  along  the 
combustor  up  to  its  end.  However,  comparing  the  heat  flux  obtained  with  those  of  cylindrical  and  dump05 
configurations  (figure  46),  we  can  see  that  it  is  of  the  same  order  of  that  of  the  first,  but  at  least  1  MW/m^  lower  than 
that  of  the  dump05  combustor. 
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Fig.  46:  Comparison  among  the  heat  fluxes  over  the  grain  surface 


As  a  consequence  also  the  regression  rate  (fig.  46)  and  the  maximum  temperature  along  the  chamber  are  lower  than 
those  of  the  DumpOS  configuration 


Fig.  47:  Comparison  among  regression  rates 
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Fig.  48:  Comparison  among  maximum  temperatures  along  the  chambers 
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Fig.  49:  Comparison  among  temperatures  along  the  chambers  (half  radius  line) 


This  effect  may  be  explained  considering  that  the  higher  dump  (1  cm)  allows  more  intense  expansion  and  thus  lower 
temperatures  (with  respect  to  the  DumpOS  configuration)  of  the  air  entering  the  chamber;  furthermore  the  Dumpi 
configuration  has  a  volume  higher  than  that  of  DumpOS.  Thus,  despite  the  better  mixing  obtainable  by  increasing  the 
dump  ratio  from  Dr  =1.1  to  Dr  =  1.2  (figures  SO,  S3)  we  cannot  reach  the  high  temperatures  obtained  in  the  others 
configurations. 
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Fig.  50:  Comparison  among  contours  of  mixing  time 


Fig.  51:  Comparison  among  mixing  times  along  the  radius  at  x  =  65  (end  of  the  grain) 
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Fig.  52:  Comparison  among  mixing  times  along  the  radius  at  x  =  85 
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Fig.  53:  Comparison  among  contours  of  mixing  time  (Dump05  vs  Dumpl) 


The  effect  of  the  increased  dump  on  the  combustion  in  the  recirculation  region  at  the  begin  of  the  grain  is  shown  in 
figure  54. 
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Fig.  54:  Comparison  among  contours  of  mixing  time  (DumpOS  vs  Dumpl-  Detail) 


As  we  can  see  a  shorter  dump  allows  shorter  mixing  times  in  the  corner,  favouring  the  flame  holding  (figures  55  and 
V56)  and  an  higher  mass  flow  rate  of  LiH  injected  in  the  reacting  layer.  That  so  enhances  the  combustion  downstream 
of  the  zone  of  the  boundary  layer  -  expansion  wave  interaction. 
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Fig.  55:  Comparison  among  contours  of  H2/O2  reaction  rate  (Dump05  vs  Dumpl-  Detail) 
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Fig.  56:  Comparison  among  contours  of  Li/02  reaction  rate  (DumpOS  vs  Dumpl-  Detail) 


The  combined  effect  of  enhanced  mixing  times  and  decreased  air  flow  temperature  on  the  rate  of  reactions  is  shown  in 
figures  57  and  58. 
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Fig.  57:  Comparison  among  contours  of  H2/O2  reaction  rate 
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In  both  cases  (Dump05  and  Dumpl  combustors)  the  H2/O2  reaction  is  mainly  developed  over  the  grain,  being  more 
intense  in  the  Dump05  configuration  and,  conversely,  more  extended  downstream  the  grain  for  the  Dumpl 
configuration. 

Analogous  considerations  holds  for  the  Li/02  reaction,  mainly  present  downstream  the  grain  and  reaching  the  core  flow 
before  the  exit  section. 


Fig.  58:  Comparison  among  contours  of  Li/02  reaction  rate 


To  complete  the  description  of  the  flow  path  in  the  Dumpl  combustor,  the  contours  of  pressure,  Mach  number  and 
mixing  time  of  both  reacting  and  non  reacting  cases  are  also  reported  in  the  following  figures  59,  60,  61. 
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Fig.  59:  Contours  of  static  pressure  (reacting  and  non  reacting  case) 
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Fig.  60:  Contours  of  Mach  number  (reacting  and  non  reacting  case) 
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Fig.  61:  Contours  of  mixing  time  (reacting  and  non  reacting  case) 
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IV  Conclusions 

The  analysis  performed  has  shown  that  LiH,  when  used  as  fuel,  offer  performance  comparable  with  those  of 
commonly  used  fuels  such  as  LCH4  and  LH2. 

During  this  preliminary  phase  we  have  in  particular  explored  their  ability  to  release  hydrogen  by  thermal 
decomposition,  showing  that  LiH  behaves  not  only  as  safe  and  compact  hydrogen  carrier  system,  but  also  as  powerful 
bi-fuels  systems,  especially  when  Li  is  injected  separately  into  the  hot  stream  of  H2/02  exhaust. 

To  overcome  the  scarcity  (or  non-existence)  of  data  concerning  LiH  combustion  at  high  temperature,  we  have 
investigated  the  behaviour  of  LiH  at  conditions  similar  to  those  expected  in  a  scramjet  combustion  chamber  by 
assuming  chemical  equilibrium.  Results  have  been  used  to  identify  main  species  and  reactions  involved  in  the  process; 
in  particular  we  have  highlighted  the  presence  of  large  amount  of  hydrogen  and  lithium  produced  by  thermal 
decomposition.  This  feature  was  stressed  comparing  the  few  theoretical  and  experimental  data  available  in  the  open 
literature  with  those  obtained  from  our  simulation;  as  a  result  we  have  demonstrated  that  at  high  temperature  the  main 
responsible  for  the  performance  of  LiH  (together  with  the  gaseous  hydrogen  produced)  is  the  liquid  lithium. 

Still  to  overcome  the  non-existence  of  data  concerning  the  Li  combustion  in  presence  of  air,  we  have 
performed  a  similar  parametric  analysis  based  on  simulations  at  chemical  equilibrium.  Results  show  that  at  temperature 
higher  than  1100  K  Li  does  not  react  with  air  nitrogen;  moreover  burning  in  oxygen  it  produces  mainly  gaseous  Li20, 
showing  at  the  same  time  many  similarities,  in  terms  of  equilibrium  composition,  with  H2/O2  combustion.  By  means  of 
a  thermochemical  analysis,  we  have  thus  compared  the  reactions  potentially  involved  in  Li  combustion,  concluding  that 
as  a  first  approximation,  the  reaction  scheme  can  be  simplified  assuming  an  one-step  global  reaction.  Furthermore,  we 
have  assumed  that  the  rate  of  this  reaction  may  be  reasonably  described  by  means  of  the  reaction  rate  proposed  by  Plane 
et  Rajasekhar.  for  the  reaction  forming  Li  superoxide. 

We  have  then  defined  the  physical  model  describing  step  by  step  the  fuel  injection  from  the  grain  surface;  we 
have  demonstrated,  by  means  of  an  order  of  magnitude  analysis,  that  almost  all  liquid  lithium  droplets  produced  by  LiH 
thermal  decomposition  (and  dragged  away  by  the  turbulent  stream),  vaporize  and  burn  before  exiting  the  combustor, 
contributing  in  a  significant  manner  to  heating  the  flow. 

To  complete  our  investigation,  we  have  also  performed  a  CFD  simulation  using  as  test  cases  three  combustor 
geometries  similar  to  those  investigated  by  Jakimovicz  et  al.  and  Ben-Arosh  et  al.  To  describe  the  fuel  injection  process 
from  the  decomposing  LiH  grain  we  have  defined  and  implemented  a  numerical  model  based  on  formulating  mass, 
momentum,  energy  and  species  sources  located  near  the  surface. 

Despite  the  limitations  due  to  the  adoption  of  a  standard  k-e  model,  and  to  the  assumption  made  about  the 
Li/02  reaction  rate,  results  are  intriguing.  The  blowing  model  adopted  describes  fuel  injection  in  a  plausible  way,  if 
not  correctly  (as  indicated  when  comparing  the  regression  rates  obtained  with  those  calculated  by  Jarymowycz  et  al.and 
Ben-Arosh).  An  intense  flame  zone  is  predicted  to  be  present  over  the  decomposing  surface  and  downstream  of  the 
grain  in  the  cylindrical  configuration;  the  flame  does  not  extinguish  and  high  temperatures  (of  order  2900  K)  are 
obtainable.  Introducing  a  sudden  cross  section  enlargement  (a  ‘dump’)  downstream  of  the  inlet  section,  mixing  in  the 
combustor  is  enhanced,  also  increasing  the  heat  released  to  the  core  flow.  However,  an  excessive  dump  ratio  decreases 
the  core  flow  temperature  (due  to  expansion  waves)  and  increases  the  combustor  volume,  lowering  the  combustor 
performance. 


56 


Grant  Number  FA8655-10-1-3091 


To  conclude  we  can  state  that  LiH  is  an  ideal  candidate  for  solid  fuelled  scramjet  applications,  behaving  as  an  high 
energy  density  bi-fuel  system  but  also  as  a  safe  and  compact  hydrogen  carrier. 
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Appendix 
A1  Evaporation 

The  simplified  evaporation  model  adopted  assumes  that  the  droplet  surface  temperature  T^,  after  a  transient  period  in 
witch  the  droplet  rapidly  heats  up  to  a  steady  temperature,  is  near  the  liquid  lithium  boiling  point  (  Tb  =  1615  K  );  this 
means  that  thanks  to  its  high  thermal  conductivity  k,  liquid  lithium  opposes  a  negligible  thermal  inertia  and,  hence,  the 
evaporation  rate  is  controlled  by  the  heat  transfer  rate  from  the  ambient  to  the  droplet  surface,  see  figure  A1 


Fig.  Al:  Droplet  and  environment  temperature  profile  adopted  in  the  model 

This  assumption  is  justified  considering  that  liquid  lithium  posses  a  high  thermal  diffusivity  (D*~5*10'^  mVs):  the 
characteristic  time  required  to  heat  a  10pm  diameter  droplet  is  thus  of  order  10'^  s,  that  is  three  order  of  magnitude  lower 
than  the  convective  time. 

Physically,  heat  from  the  environment  supplies  the  energy  necessary  to  vaporize  the  liquid  fuel,  and  the  fuel  vapor  then 
diffuses  from  the  droplet  surface  into  the  core  gas.  The  mass  loss  causes  the  droplet  radius  to  shrink  with  time  until 
the  droplet  is  completely  evaporated.  This  is  a  good  approximation  especially  for  the  scramjet  combustion  chamber 
where  the  expected  maximum  temperature  is  very  high  (Tc  about  3500  K);  thus  the  problem  is  reduced  to  determining 
the  mass  flowrate  of  the  fuel  vapor  from  the  surface  at  any  instant  in  time.  Knowledge  of  this  will  then  enable  us  to 
calculate  the  droplet  radius  as  a  function  of  time,  and  droplet  lifetime. 

In  summary,  the  model  include  the  following  assumptions,  commonly  adopted  in  evaporation  problems  because  they 
lead  to  great  simplification  while  maintaining  a  reasonable  agreement  with  experimental  results;’* 
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-  The  droplet  evaporates  in  a  quiescent,  infinite  medium. 

-  The  evaporation  process  is  quasi-steady.  This  means  that  at  any  instant  in  time  the  process  can  be  described 

as  if  it  were  in  steady  state.  This  assumption  eliminates  the  need  to  deal  with  partial  differential  equations. 

-  The  fuel  is  a  single-component  liquid  with  zero  solubility  for  other  gases  present. 

-  After  a  short  transient,  the  droplet  temperature  is  uniform,  and,  furthermore,  the  temperature  is 
assumed  to  be  the  lithium  boiling  point,  Tj  =  Tb.  In  many  problems,  the  transient  heating  of  the  liquid 
does  not  greatly  affect  the  droplet  lifetime,  and  more  rigorous  calculations  show  that  the  droplet  surface 
temperature  is  only  slightly  less  than  the  liquid  boiling  point  in  combustion  environments.  This  assumption 
eliminates  the  need  to  solve  a  liquid-phase  (droplet)  energy  equation,  and,  more  importantly,  eliminates  the 
need  to  solve  the  fuel  vapor  transport  (species)  equation  in  the  gas  phase.  Implicit  in  this  assumption  is  that 

Too  >  Tb. 

-  We  assume  binary  diffusion  with  a  unity  Lewis  number  (a  =  D).  This  permits  us  to  use  the  simple  Shvab  - 
Zeldovich  treatment  for  the  energy  equation. 

-  We  also  assume  that  all  thermophysical  properties,  such  as  thermal  conductivity,  density,  and  specific  heat, 
are  constant.  Although  these  properties  may  vary  greatly  as  we  move  through  the  gas  phase  from  the  droplet 
surface  to  the  faraway  surroundings,  constant  properties  allow  a  simple  closed-form  solution.  In  the  final 
analysis,  a  judicious  choice  of  mean  values  allows  reasonably  accurate  predictions  to  be  made. 

Using  these  we  can  find  the  vaporized  lithium  mass  flowrate  focusing  only  on  the  gas  phase  analysis;  thus  we  can  write 
the  mass  conservation  equation  as: 


m  =  p  v^4  71  =  constant 


(Al.l) 


where  ih ,  r  ,  n^  are  respectively  the  gaseous  lithium  mass  flowrate  (constant  because  of  steady  state  evaporation 
assumptions),  its  density  and  the  bulk  flow  velocity.  Thus  the  conservation  of  energy  is: 


dr 


dr 


47ik  dr 


(A1.2) 


where  we  have  assumed  that  no  reaction  are  present.  This  equation  requires  the  following  boundary  conditions: 
far  from  the  droplet :  T  (r  ^  oo)  =  Tstream 
at  the  drop  surface:  T  =  Tboii 
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Finally  we  can  write  the  droplet  -  gas  phase  interface  energy  balance  as 


Q, 


cond 


A  ,  2cIT 
“''k.r.- 


m(h 


vap 


K)  +mAh 


heating 


(A1.3) 


were  m(h^^p- h,;J  is  the  lithium  latent  heat  of  vaporization  and  is  the  heat  required  to  reach  the 

lithium  boiling  temperature. 

Solving  equation  (A1.3)  we  obtain  the  temperature  distribution  as  a  function  of  the  droplet  diameter;  thus  differentiating 
this  relation,  substituting  in  (A1.3)  and  solving  for  lil  we  obtain; 


4  71  k„  r 


m : 


g  s 


ln(B,+l) 


(A1.4) 


were  Bq,  sometimes  referred  to  as  the  Spalding  number  or  transfer  number,  is: 


Cpg  (T„  Tj,oii) 


_  pg _ 


(A1.5) 


The  droplet  lifetime  can  thus  determined  by  following  its  radius  history  by  writing  an  equation  stating  that  the  rate  of 
the  droplet  mass  consumption  equals  the  vaporized  mass  flowrate: 


dt 


=  -  rh 


(A1.6) 


and,  since  =  PlV  =  71  D^/6  we  can  write  a  relation  that  is  commonly  expressed  in  terms  of  rather  than  D; 


dP" 

dt 


8  k 


Pl  c 


L  ''pg 


ln(B,+l) 


(A1.7) 


and,  solving; 
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D"(t)  =Do  -Kt 


(A1.8) 


8k 


U  JV  /  X 

with  K  =  -  In  +lj  defined  as  the  evaporation  constant  (figure  A1.2)  and  Dq  the  initial  droplet  diameter. 


Pl  c 


L  ''Pg 


Fig.  A1.2:  law 


We  can  find  the  time  it  takes  a  droplet  of  given  initial  size  to  completely  evaporate,  i.e.,  the  droplet  lifetime  t^,  by 
setting  D^(t(])  =  0; 


t 


d 


(A1.9) 


For  the  case  under  examination  this  characteristic  time  was  calculated  as  a  function  of  the  temperature  with  the  drop 
diameter  as  parameter;  this  choice  was  made  because  during  the  residence  time  in  the  combustion  chamber  the  droplet 
are  exposed  to  a  temperature  field,  depending  on  where  along  the  chamber  axis  they  are  released  and  heated.  Due  to  the 
temperature  dependence  of  the  air  thermal  conductivity  kg  and  constant  pressure  specific  heat  Cpg,  an  average  value 
was  defined  for  every  ambient  temperature;^*’’^ 


Tm  -  (Tboil  +  T*,)  /  2 


(ALIO) 


Thus  kg  and  Cpg  were  calculate  as  average  values  of  the  air  and  lithium  properties  calculated  at 


^Pg  “  ^p(Li)(Tin) 


(Al.ll) 


and 
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kg=0.4k(u)(TJ  +  0.6k,,(TJ 


(A1.12) 


Lithium  and  air  properties  were  calculated  using  the  polynomials  in  the  CEA600  database;  results  are  shown  in  figure 
A1.3 


- 1  micron 

—  10  micron 

- 50  micron 

- Conv.  Time 


Flow  Temperature  [K] 


Fig.  A1.3:  Evaporation  rate 

As  we  can  see  the  lifetime  of  droplets  having  diameters  of  few  microns  is  of  the  order  of  ten  millisecond  at  low 
temperatures  (typical  of  the  initial  transient  phase  or  far  from  the  reacting  zone);  however  it  is  important  to  highlight 
that  these  values  are  relatives  to  a  quiescent  and  non  reacting  environment  so  neglecting  the  convective  heat  transfer 
and  the  coupled  effects  of  the  fuel  mass  driven  diffusion  and  combustion.  These  will  be  analyzed  in  the  next  section. 


Al.l  Droplets  evaporation  and  burning’^ 

We  can  extend  the  model  developed  in  the  preceding  section  including  a  spherically  symmetric  diffusion  flame  that 
surrounds  the  droplet.  While  retaining  the  assumption  of  a  quiescent  environment  and  spherical  symmetry  and  short 
heat-up  time  made  in  the  initial  development,  we  now  couple  in  the  analysis  the  mass  and  heat  transfer.  The  method  is 
thus  based  on  the  following  assumptions: 

-  The  burning  droplet,  surrounded  by  a  spherically  symmetric  flame,  exists  in  a  quiescent,  infinite  medium.  There 
are  no  interactions  with  any  other  droplets,  and  the  effects  of  convection  are  ignored.  This  reasonable,  as  the 
droplets  are  stripped  from  the  surface,  so  that  the  relative  velocity  is  either  zero  or  very  small  compared  to  the 
core  gas  velocity. 
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-  As  in  our  previous  analysis,  the  burning  process  is  quasi-steady. 

-  The  fuel  is  a  single-component  liquid  with  zero  solubility  for  gases.  Phase  equilibrium  prevails  at  the  liquid- 
vapor  interface. 

-  The  pressure  is  uniform  and  constant. 

-  The  gas  phase  consists  of  only  three  "species:"  lithium  vapor,  air  and  combustion  products.  The  gas-phase 
region  is  divided  into  two  zones.  The  inner  zone  between  the  droplet  surface  and  the  flame  contains  only 
lithium  vapor  and  products,  while  the  outer  zone  consists  of  air  and  products.  Thus,  binary  diffusion  prevails 
in  each  region. 

-  Lithium  and  air  react  in  stoichiometric  proportions  at  the  flame.  Chemical  kinetics  is  assumed  to  be  infinitely 

fast,  resulting  in  the  flame  being  represented  as  an  infinitesimally  thin  sheet. 

-  The  Lewis  number  is  unity. 

-  Radiation  heat  transfer  is  negligible. 

-  The  gas-phase  thermal  conductivity,  kg,  specific  heat,  Cpg  and  the  product  of  the  density  and  mass  diffusivity, 
pD  are  all  constants. 

-  The  liquid  fuel  droplet  is  the  only  condensed  phase. 


A  schematic  representation  of  the  model  is  shown  in  figure  A  1.4  where  we  have  outlined  the  process  driving  factors 
(that  is  the  temperature  and  mass  fraction  gradients)  as  a  function  of  the  radial  distance  from  the  droplet  core. 
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Fig.  A1.4:  Temperature  and  mass  fractions  as  a  function  of  the  radial  distance 


This  model  thus  involves  the  simultaneous  solution  of  the  following  set  of  equation: 


-  An  equation  governing  the  lithium  conservation  in  the  inner  region  (that  is  rs  <  r  <  rf  with  rs  and  rf  respectively 
the  surface  and  flame  sheet  radial  coordinates) 


-  An  analogous  equation  governing  the  air  species  conservation  in  the  outer  region  (rf  <  r  <  r^) 


-  The  energy  conservation  at  the  droplet  liquid-vapor  interface 


-  The  energy  balance  at  the  flame  sheet 


-  The  Clausius  -  Clapeyron  equation  governing  the  liquid-vapor  equilibrium  at  the  interface. 


Solving  the  corresponding  equation  system  we  obtain  the  vaporized  lithium  mass  flowrate: 


m= - ^ — ln(B„^+l)  (A1.13) 

^pg 
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were  Ahj.  is  heat  of  combustion  and  v  is  the  lithium-air  stoichiometric  mass  ratio;  thus  the  transfer  number  Bo,q  is 
defined  as: 


Ah„ 


+  -TJ 


B  =  - 

Ah,„  -f  Ah 


heating 


(A1.14) 


the  burning  rate  constant  is: 


8  k,  I  . 

K= - ^ln(B  +l)  (A1.15) 

Pl  Sg 

and  the  droplet  lifetime: 

td  =  (A1.16) 

Is. 

As  with  the  pure  evaporation  problem,  the  appropriate  values  of  the  gas  thermal  conductivity  kg  and  constant  pressure 
specific  heat  Cpg  were  calculated  defining  an  average  temperature  T,i;  since  Tf  =  3300  K,  the  flame  temperature  of  the 
lithium-air  combustion  calculated  using  the  NASA  CEA2  code,  T,j  is  defined  as: 

T„,=  (Tf+T,)/2  (A1.17) 

with  Ts  the  droplet  surface  temperature.  Results  are  shown  in  figure  A  1.5  where  the  droplet  lifetime  as  a  function  of 
the  environment  temperature  and  the  droplet  diameter  was  compared  with  the  convective  time. 


10  -| 

(A 

£ 

1  - 

o 

£ 

“S) 

0,1  - 

-1 

Q. 

p 

0,01  - 

D 

0,001  -1 

-1  micron 
-10  micron 
-50  micron 
-Conv.  Time 


1800  2200  2600  3000  3400 

Chamber  Temperature 


3800 


Fig.  A1.5:  Droplets  lifetime 
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A  more  accurate  evaluation  of  the  gas  characteristic  was  obtained  using  values  of  kg  and  Cpg  calculated  using  CEA600; 
results  are  shown  in  figure  A1.6 
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Fig.  A1.6:  Droplets  lifetime  (data  obtained  with  CEA600) 


A1.2  Extension  to  a  convective  environment 

The  last  step  of  our  investigation  on  lithium  liquid  droplets  evaporation  and  combustion  is  to  incorporate  the 
convection  effect  into  the  model  developed.  We  follow  the  approach  suggested  by  Turns,’*  based  on  the  “film  theory”. 
The  difference  between  this  last  and  the  preceding  model  consists  in  the  replacement  of  the  boundary  conditions;  that 
is  the  same  boundary  conditions  used  at  infinity  are  now  moved  inward  to  a  distance  5m  for  species  and  5t  for  energy, 
see  figure  A  1.7. 
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Fig  A1.7:  Effect  of  the  convective  environment 

In  figure  A1.7,  the  film  radius  steepens  the  concentration  and  temperature  gradients,  so  enhancing  mass  and  heat 
transfer  rates  at  the  droplet  surface.  This  means  that  convection  increases  the  drop  burning  rate  and  lowers  the 
characteristic  time  with  respect  to  the  preceding  case  (quiescent  environment).  The  film  radii  are  defined  as: 


^  _  Nu 
~  Nu-2 


(A1.18) 


were  the  Nusselt  number  physically  represents  the  dimensionless  temperature  gradient  at  the  droplet  surface  and 


5m  Sh 
r,  Sh  -  2 


(A1.19) 


were  the  Sherwood  number  Sh  is  the  dimensionless  concentration  at  the  surface.  We  can  observe  that  for  a  quiescent 
environment  Nu  =  2;  in  the  limit  5t  ^  co,  results  agree  with  the  boundary  conditions  used  into  the  preceding  case. 
Assuming  Sh  =  Nu  (remembering  that  the  Lewis  number  Le  =  1)  a  correlation  to  evaluate  Nu  is: 


Nu  =  2  + 


0.555 

1  +  1.232/(RePr"^^)]‘^" 


(A  1.20) 


70 


Grant  Number  FA8655-10-1-3091 


were  Re  is  based  on  the  droplet  diameter  and  the  relative  velocity,  while  the  thermophysical  properties  can  be 
calculated  as  in  (ALIO  ;  A1.12). 

Finally  the  burning  rate  constant  (A1.15)  that  can  be  extracted  from  this  theory  is: 


K  =  Nu 


4  k, 

Pl  Sb 


(A1.21) 


The  impact  of  these  more  realistic  assumption  (considering  the  high  flow  speed  into  the  chamber)  can  be  shown  by 
examining  what  happens  to  a  droplet  having  a  larger  diameter  than  allowed  by  the  analysis  developed  in  the  preceding 
sections.  Assuming  in  fact  the  existence  of  100  )J,m  diameter  droplets,  we  can  evaluate  the  time  employed  by  one  of 
these  drops  to  reach,  starting  from  rest,  a  velocity  corresponding  to  60%  of  the  flow  speed:  this  time,  called  relaxation 
time  Xp  accounting  for  the  stream  turbulence,  is  defined  by  modifying  the  laminar  relaxation  time 
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with  the  term  in  squared  brackets,  becoming  thus:’® 
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where  Rej,  Pu  ,  D  are  respectively  the  Reynolds  number  of  the  droplet,  the  liquid  lithium  density  and  the  droplet 
diameter.  Using  the  same  numerical  values  in  the  preceding  section,  the  droplet  relaxation  time  is  of  the  order  of  10'"^ 
seconds,  that  is  smaller  than  the  convective  time;  incidentally  we  can  observe  that  calculated  the  corresponding 
acceleration,  the  distance  covered  by  the  drop  during  a  time  interval  equal  to  the  convective  time  is  about  50  cm, 
showing  that  there  is  enough  time  for  heating  and  combustion.  Now  if  we  calculate  with  (A  1.21)  and  (A  1.23)  the  Nu 
and  K  at  the  initial  instant  to  =  0  and  after  ti  =  Xp  seconds,  we  obtain  for  the  droplet  lifetime  the  results  presented  in 
figure  A  1.8: 
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Fig  A1.8:  100  |im  droplet  lifetime  in  convective  environment 

Thus,  despite  the  droplet  size,  its  lifetime  at  the  initial  instants  of  the  motion  is  of  the  same  order  of  the  relaxation 
time;  and  after  a  time  interval  corresponding  to  the  relaxation  time,  the  drop  lifetime  is  still  one  half  of  the 
characteristic  convective  time.  The  lifetime  is  even  shorter  obviously  when  considering  droplet  having  diameter 
smaller  that  100  pm. 
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